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Abstract. We combine single- and two-plroton interference procedures for 
characterizing any multi-port linear optical interferometer accurately and precisely. 
Accuracy is achieved by estimating and correcting systematic errors that arise due to 
spatiotemporal and polarization mode mismatch. Enhanced accuracy and precision are 
attained by fitting experimental coincidence data to curve simulated using measured 
source spectra. We employ bootstrapping statistics to quantify the resultant degree of 
precision. A scattershot approach is devised to effect a reduction in the experimental 
time required to characterize the interferometer. The efficacy of our characterization 
procedure is verified by numerical simulations. 
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Linear optics is important in quantum computation and communication. The simulation 
of a linear optical interferometer is computationally hard classically subject to reasonable 
conjectures [1]. Single-photon detectors and linear optical interferometers allow 
for efficient universal quantum computation via linear optical quantum computing 
(LOQC) [2]. Linear optics can simulate the quantum quincunx [3] and quantum random 
walks [4]. Linear optics coupled with laser-manipulated atomic ensembles enables long¬ 
distance quantum communication [5]. A wide class of communication protocols can be 
realized with coherent states and linear optics [6]. 

Recent advances in photonic technology including photonic circuits on silicon chips [7- 
11], noise-free high-efficiency photon number-resolving detectors [12-16], high-fidelity 
single-photon sources [17-20] have engendered the experimental implementation of multi- 
port linear optical interferometry. Reconfigurable linear optical interferometers that can 
perform arbitrary unitary transformations have been demonstrated [21-23]. 

The accurate and precise characterization of linear optics is important in quantum 
information processing tasks such as BosonSampling, LOQC and quantum walks. 
BosonSampling involves sampling from the output photon coincidence distribution of an 
interferometer on single photon inputs to each mode. Sampling from this distribution is 
computationally hard classically but is easy with a linear-optical interferometer. The 
classical hardness of the BosonSampling problem crucially depends on the error in the 
linear optical interferometer [24]. Similarly, the practical applications of BosonSampling, 
in quantum metrology and in the computation of molecular vibronic spectra, rely on the 
accurate implementation and characterization of linear optics [25,26]. 

Accurate and precise characterization is important in LOQC because a high success 
probability of the employed non-deterministic linear-optical gates relies on implementing 
the desired gates with high fidelity [27]. Furthermore, linear interferometers used 
in photonic quantum walks, which display strong non-classical correlations, require 
accurate characterization especially if quantum walks are employed for solving classically 
hard problems [28-30]. In other words, that accurate and precise characterization of 
interferometers enables a verifiable quantum speedup of linear-optical protocols over 
classical computers. 

Classical-light procedures [31,32] for linear optics characterization are unsuitable 
for Fock-state based experiments because the interferometer parameters change during 
the coupling and decoupling of classical light sources and of homodyne detectors at the 
interferometer ports. This change could result from drift of interferometer parameters in 
the time required to couple sources and detectors or as the result of mechanical process of 
coupling itself. Characterization procedures that rely on Fock-state (rather than classical- 
light) inputs are thus more desirable in BosonSampling and LOQC implementations; such 
procedures would enable interferometer characterization without altering the experimental 
setup and would thus be accurate. 

The Laing-O’Brien procedure [33] uses one- and two-phot.ons for characterizing 
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linear optical interferometers and is stable to the length scale of a photon packet. This 
procedure assumes perfect matching in source field and large-number statistics on the 
detected photons. Hence, implementations of this procedure are inaccurate due to 
spatiotemporal and polarization mode mismatch in the source field and imprecise due to 
shot noise. 

We aim to devise an accurate and precise procedure that uses one- and two-photons 
for the characterization of linear optical interferometers and to devise a rigorous method 
to estimate the standard deviation in the linear optical interferometer parameters [34], 
Furthermore, we aim to provide a correct alternative to the x 2 -tesf, which has been 
used to estimate the confidence in the characterized interferometer parameters in current 
BosonSampling implementations [11, 35, 36]f. 

Here we devise a procedure to characterize a linear optical interferometer accurately 
and precisely using one- and two-phot.on interference. Four strengths of our approach 
over the Laing-O’Brien procedure [33] are that our procedure (i) accounts for and 
corrects systematic error from spatial and polarization source-field mode mismatch via 
a calibration procedure (Section 3); (ii) increases accuracy and precision by fitting 
experimental coincidence data to curve simulated using measured source spectra 
(Section 3); (iii) accurately estimates the error bars on the characterized interferometer 
parameters via a bootstrapping procedure (Section 4); and (iv) reduces the experimental 
time required to characterize interferometers using a scattershot procedure (Section 5). 

2. Background 

This section provides the background for our one- and two-photon characterization 
procedure. The action of a multi-port linear optical interferometer on single photons 
entering one or two input ports and vacuum entering the other ports is detailed. 
Specifically, we calculate the probability of detecting a photon at a given output port when 
a single photon is incident at a given input port. The section concludes with expressions 
for the probability of detecting a coincidence measurement when two controllably delayed 
photons are incident on the interferometer. 

2.1. Action of a linear optical interferometer 

In this subsection, we define linear optical interferometers by their action on single 
photons. We parameterize the unitary transformation effected by an interferometer and 

f The % 2 -test [37-39] is used to quantify the goodness of fit between probability distribution functions 
of two categorical variables, which can take a fixed number of values. Coincidence-count curves and 
visibilities are not probability distribution functions of categorical variables, but rather are collections 
of many categorical variables (variables that can take on one of a fixed finite number of possible values), 
one variable corresponding to each time-delay value chosen in the experiment. Hence, quantifying 
the goodness of fit between two coincidence curves using the y 2 -test is incorrect. This incorrectness 
undermines the claim that the data are consistent with quantum predictions and disagree with classical 
theory [11,36] and leaves the choice of unitary matrices [35] unjustified. 
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present our treatment of losses and dephasing at the interferometer ports. 

Consider a single photon entering the i-th mode of an m-mode interferometer. The 
monochromatic photonic creation and annihilation operators acting on the i-th and the 
j'-th ports obey the canonical commutation relation§ 

a.i(uji),a](uj2) 

for positive real frequencies uq ,uu 2 . The state of a single photon entering the i-th mode is 

/ OO 

dufi(u)a\(uj) | 0 ), ( 2 ) 

-OO 

where fi(co) is the normalized square integrable spectral function, |0) is the m-mode 
vacuum state. The state of two photons entering modes i and j i of the interferometer 

is 

/ oo roo 

dcai / dca 2 /* (uq) fj (u> 2 )aj(oq) a j (ca 2 ) 10) (3) 

oo J — oo 

with exchange symmetry holding if fjfuj) = /j(ca). One- and two-photon states are 
transformed into superpositions of one- and of two-photon states respectively under the 
action of the linear interferometer. 

We treat linear interferometers as unitary quantum channels acting on the state of 
the incoming light. The interferometer transforms the photonic creation and annihilation 
operators according to 

m 

a ) M Vi i M ( 4 ) 

and its complex conjugate, where V{co) is the transformation matrix of the interferometer. 
Photon-number conservation imposes unitarity 

V*{u)V{u) = 1 (5) 

of the transformation matrix V{co) for all real u. In general, the elements {Vijico)} of 
the transformation matrix depend on the frequency of transmitted light. We assume 
that the spectral functions of the incoming light are narrow compared to frequencies over 
which the entries {W,-} change noticeably and thus treat V to be frequency-independent. 

If only Fock states are incident at the interferometer and only photon-number- 
counting detection is performed on the outgoing light, then the measurement outcomes 
are invariant under phase shifts at each input and output port. That is, interferometer 
V = DiVD\ produces the same measurement outcome as V for any diagonal unitary 
matrices Di and D 2 . Mathematically, if Di,D 2 are diagonal unitary matrices, then 

V~v <=> V = D 1 VD\ ( 6 ) 

is an equivalence relation. Members of the same equivalence class defined by this 
equivalence relation produce the same number-counting measurement outcomes on 
receiving Fock-state inputs. 

§ Two monochromatic photons are distinguishable based on the ports that they occupy and on their 
respective frequencies uq and w 2 . 


— 8ij5(u\ — CP 2 )1 


( 1 ) 
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Figure 1: Schematic diagram of the interferometer. U effects a unitary transformation on 
a multimode state of light. The dotted lines represent the couplings of the interferometer 
with light sources and detectors. The beam splitters at the input and output modes 
model the linear losses because of imperfect coupling and detector inefficiency. The 
vacuum input to these beam splitters is not shown. One of the beam splitter outputs 
enters the interferometer while the other one is lost. The triangles represent the random 
dephasing at the input and output ports. The dashed box labelled U lossy represents the 
combined effect of the dephasing, the losses and the unitary interferometer. 


Each equivalence class can be represented by a unique matrix U whose first row 
and first column consist of real elements. The complex matrix entries of the class 
representative U ~ V are 

{Uij = t t] e' (kj : tij e R + , 6 i:j e (-7r, 7r], 0 a = 0, 9 Xj = 0 Vz, j e {1,2,..., m }} . (7) 

The constraints 9n = 0, = 0 Vz G {1,2,..., m} on the input and output phases of the 

transformation matrix are obeyed in the following parameterization of U 

U = L x A x M, 
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Thus, the values {A*}, {a^}, {%}, {p,j} completely parametrize the class representative 
matrix U. 

The input and output ports of the interferometer are amenable to time-dependent 
linear loss and dephasing. We model losses using parameters zq- and zq, which are 
the respective probabilities of transmission at the input mode j and output mode z. 
Dephasing is modelled using parameters and 0,, which are the arbitrary multiplicative 
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phases at the input and output ports. Hence, the actual transformation effected by the 
interferometer is given by the matrix U lossy , which has matrix elements 

jj\o Ssy = 

= e 1 ^ -s/Ki^/Xi aije l9i i yfJTjy/ife 1 ^. (9) 

Figure 1 depicts the relation between the representative matrix U and the actual 
transformation £/]° hSy that is effected by the interferometer. 

This completes our parameterization of the linear optical interferometer. Our 
characterization procedure employs one- and two-photon inputs to estimate the values of 
parameters {A*}, {a^}, {%}, {p-j} of (9). In the next subsection, we recall the expectation 
values of measurements performed on interferometer outputs when one- and two-photon 
states are incident at the input ports. 

2.2. One- and two-photon inputs to linear optical interferometer 



Figure 2: Schematic diagram of single-photon counting at the output of an interferometer 
when single photons are incident at one input port. The star symbol represents a source 
of single-photon pairs. Single photons are incident at one of the input ports while vacuum 
state is input to the remaining input ports (not shown in figure). The semicircles at 
the output ports represent single-photon detectors and the circles with the included ff 
represent the photon-counting logic connected to the detectors. 


Our characterization procedure employs single-photon counting to estimate the 
amplitudes {a^} of the representative matrix U entries. The arguments {%} of 
U are estimated using two-photon coincidence counts. In this subsection, we give 
expressions for one- and two-photon transmission probabilities, which are employed in 
our characterization procedure (Section 3). 











































Accurate and precise characterization of linear optical interferometers 


7 


We first consider the case of single-photon transmission. The interferometer 
transforms the single-photon input state ( 2 ) to the state at the output ports according 
to (9). A photon is detected at the i-th output port with a probability 


P = 

1 13 


rrlossy 


\ 2 

— ij flj Vj. 


( 10 ) 


when a single-photon is incident on the j-tli input port. 

Whereas the values of {ay,} are estimated using single photon counting, {%} values 
are estimated using two-photon coincidence measurement. We now present probabilities 
of detecting two-phot.on coincidence at the interferometer outputs when cont.rollably 
delayed pairs of photon are incident at the input ports. 

If a controllably delayed photon pair is incident at input ports j and f, then the 
probability Cn'jj>(r) of coincidence measurement at detectors placed at output ports i 
and i' is 


Cu'jji ( t ) = KiKi'VjVji (tijtlj, + j duiduilfjMfj'MY 

+ 2 tijtijdyjti'j' J daydtc 2 fj(ui)fj' (cu 2 ) fj (^2 )fj / ( UJ i) 

X COS (l0 2 t — 0J\T + 6ij — dijI — 9i'j + Oi'jf) 

On substituting according to (7), we obtain [40] 


( 11 ) 


C ii' jj '( T ) =KiKi>W'H j Hi'l / j l' j ' { a ij a i'j' +a2 if a h) j dcci dcca \fj(uJi)ff {u 2 ) I s 
+ 27 OLijOiij'OLi'jOti'j' J duidu 2 fj{u 1 )f jl (u 2 )fj{u} 2 )fj'{ui) 

X COS (u 2 T — 0J\T + dij — dij / — di'j + di'j ') 


( 12 ) 


where r is the time delay between the two photons, fj(u>), fj'ipj) describe the spectrum 
of light just before it enters the detectors and 7 is the mode-matching parameter, which 
we described in the remainder of this section. 

Two-phot.on coincidence probabilities (12) depend on the mode matching in the 
source field. Spatial and polarization mode mismatch is quantified by the mode-matching 
parameter 7 [40]. Perfectly indistinguishable light sources, such as light from a single¬ 
mode fibre, have relative mode matching 7 = 1 whereas 7 = 0 indicates that the sources 
are completely distinguishable. Figure 4 depicts how imperfect mode matching, i.e., 
7 < 1, alters the observed t.wo-phot.on coincidence counts. Our calibration procedure 
estimates and accounts for imperfect mode matching, which is assumed to be constant 
over the runtime of the characterization experiment. 

The calculation of the expected coincidence probabilities as a function of the 
time delay between the photons is detailed in Algorithm 1. The next, section 
describes how single-photon transmission probabilities ( 10 ) and two-photon coincidence 
probabilities ( 12 ) are used for characterizing the linear optical interferometer. 
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Figure 3: Schematic diagram for coincidence measurement the interferometer output 
when single-photon pairs are incident on two different input ports of an interferometer. 
The star symbol represents a source of single-photon pairs and the semicircles at the 
output ports represent single-photon detectors. The coincidence logic, which is depicted 
by 0, counts two-photon coincidence events at the detectors. 



Figure 4: Plots of coincidence probability versus time delay for different values of 7 for 
a lossless balanced beam splitter. The time delay r is in units inverse special width of 
incoming photons. 


3. Characterization of linear optical interferometer 

In this section, we describe our procedure to characterize linear optical interferometers. 
The outline of this section is as follows. Subsection 3.1 describes the experimental 
data required by our characterization procedure. This experimental data are processed 
by various algorithms to determine the transformation matrix ( 8 ). The algorithm to 
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determine the amplitudes {ay,} of the transformation-matrix elements is presented in 
Subsection 3.2. In Subsection 3.3, we describe the calibration of the source held by 
determining the mode-matching parameter 7. The estimation of {%} using two-photon 
interference is detailed in Subsection 3.4. Maximum-likelihood estimation is employed 
to find the unitary matrix U that best fits the calculated {ay,}, {%} values and serves 
as the representative matrix ( 8 ). We discuss the calculation of the best-fit unitary 
representative matrix in Subsection 3.5. 

3.1. Experimental procedure and inputs to algorithms 

Our characterization procedure relies on measuring (i) the spectral function fj of the 
source light, (ii) single-photon detection counts, (iii) two-photon coincidence counts 
from a beam splitter and (iv) two-photon coincidence counts from the interferometer. 
The measurement data constitute the inputs to our algorithms, which then yield the 
representative matrix. Before presenting the algorithms, we detail the experimental 
procedure and the inputs received by the algorithm in this subsection. 

We characterize the spectral function /(to*) of the incoming light for a discrete set 
Q = {o;i,a; 2 ,... ,07,} of frequencies. The integer k of frequencies at which the spectral 
function is characterized is commonly equal to the ratio of the bandwidth to the frequency 
step of the characterization device. The characterized spectral function ficof) is used to 
calculate the coincidence probabilities as detailed in Algorithm 1. 

Algorithm 1 COINCIDENCE: Calculates the expected coincidence rate for two- 
photon interference for a given 2 x 2 submatrix of an arbitrary SU(m) transformation. 

Input: 

• k, f2 = ( 07 , 07 ,... u>k-i,uJk} £ (R + ) fc > Frequencies at which fi, / 2 are given. 

• / 1 , / 2 : O —> R + > measured spectra. 

• £, T — {77, r 2 ,...,Tf }6 (RU0) f > Time delay values. 

• A G- {cc^, tty/, a t] , } G (R + U 0 ) 4 > Amplitudes of 2 x 2 submatrix of A ( 8 ). 

• 0 G- Ojj , $iji, di'j , 9 iijt G (— 7 r, 7 r] > Phases of 2 x 2 submatrix of A ( 8 ). 

• 7 G [0,1] > Mode-matching parameter of photon source. 

Output: 

• C : T —y R + > Two-photon coincidence probabilities correct up to 

multiplicative factor. 

1: procedure Coincidence^, Q, / 1; / 2 , £, T, A, 0 ,7) 

2 : for r in T do 

3 : C(t) <- Integrate [F(A, < 3 >, /1, / 2 ,7,07, uij, r), {u a ett,uj b e 12 }]. 

> Numerically integrate RHS of (12) over 07,07 with K. t = K t > = Uj = vy = 1. 

4: end for 

5: return C 

6: end procedure 


The amplitudes {op} are determined by impinging single photons at the 
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interferometer and counting single-photon detections at the outputs. Single-photon 
counting is repeated multiple (B G Z + ) times in order to estimate the precision of the 
obtained {ay,} values. Specifically, the number 

N ijbj G Z + : i,j G {1 e (13) 

of single-photon detection events are counted at all m output ports {i} for single 
photons impinged at the j-th input ports in the bj -th repetition. The counting is then 
performed for each of the input ports j G {1,..., m } of the interferometer. Algorithm 2 
uses bj £ { 1 , • • •, -B}} values to estimate a tJ and the standard deviation of the 

estimate. The experimental setup for {ay,} measurement is depicted in Figure 2. 

Arguments {%} are calculated by fitting curves of measured coincidence counts 
to curves calculated using measured spectra according to (12). Appendix B elucidates 
the inputs and outputs of the curve-fitting procedure, such as the Levenberg-Marquardt 
algorithm [41,42], employed by our algorithms. Before calculating {%}, we calibrate the 
source field for imperfect mode matching by measuring coincidence counts on a beam 
splitter of known reflectivity. Cont-rollably delayed single-photon pairs are incident at¬ 
tire two input ports of the beam splitter and coincidence counting is performed on the 
light exiting from its two output ports. Algorithm 3 details the estimation of 7 using 
coincidence counts C cal (r) for time delay r between the incoming photons. 

The absolute values and the signs of the arguments {% G (—7r, 7 t] } are calculated 
separately. To estimate the absolute values {1 6 t j \} of the arguments, pairs of single 
photons are incident at two input ports 1 and j G {2, ...,m} and coincidence 
measurement is performed at two output ports 1 and i G {2, ...,m}. The choice 
of the input and output ports labelled by index 1 is arbitrary. The signs 

f -1 if'% < 0 , 

sgn % = <j 0 if dij — 0, (14) 

[l if 6ij > 0 

of the arguments are estimated using an additional (m — l ) 2 coincidence measurements. 
Algorithm 6 details the choice of input and output ports for estimating {sgn^j,,}. A 
schematic diagram of the experimental setup for {%} estimation is presented in Figure 3. 

3.2. Single-photon transmission counts to estimate {ay,} (Algorithm 2) 

Now we present our procedure to estimate {ay,} values using single-photon counting. 
Single-photon transmission probabilities are connected to the amplitudes {ay,} according 
to the relation P yJ = K^AjoA/ij/y, (10). Although the {ay,} values can be calculated 
from single-photon transmission counts, the factors {A;}, {/y,} cannot. The transmission 
probabilities depend on the products of the factors {A*}, {/y,} and the loss terms {«;*}, {zy,}, 
so {Aj}, {pj} cannot be measured without prior knowledge of the losses. The loss terms 
are usually unknown and can change between experiments. Hence, we calculate the 
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values of {a tl } from single-photon measurements and choose {A,} and {j-ij} such that 
U = LAM is unitary. 

The amplitudes {ay,} are determined by estimating transmission probabilities. The 
probabilities Pu, Pn, Pij, Pij of single-photon detection at output ports 1 ,i when single 
photons are incident at input ports 1 ,j are expresses in terms of the a VJ values according 
to 

P\ 1 Pij \r lAiO!ii/iiSi| \r[AjSj 

P\jP %i |ri \\Oi\jp,jSj | \riXiCXupiSi 

The probabilities Pn, Pn, Pij, Pij are estimated by counting transmitted photons. The 
definition (8) of a tJ implies that an = an = a\j = 1. Hence, the values of a l3 are 
connected to the single-photon transmission probabilities according to 


^ 11 




(15) 


PnPij 

a ‘ i ~' l Pi J Pa’ 


(16) 


which is independent of the losses at the input and the output ports. 

The transmission probabilities P tJ are estimated by counting transmitted photons 
as follows. The estimated values of {a t j} are random variables that are amenable 
to random error from under-sampling and experimental imperfections. Thus, data 
collection is repeated multiple times. For accurate estimation of cpj and its standard 
deviation daij , the number B of repetitions is chosen such that the standard deviation of 
{Nij^ : bj G {1,..., B}} converges in B for all i, j G {1,..., m}. The mean and standard 
deviation of {N t] b :i '■ bj G {1,..., B}} converge for large enough B if the cumulants of 
the distribution are finite [43]. 


Algorithm 2 AmplitudeEstimation: Uses single-photon detection counts to 
calculate the amplitudes of the complex entries of the transformation matrix. • 
represents our estimate of •. 

Input: 

• m G Z + , > Number of modes of interferometer. 

• N ijbj : {1,..., m} x {1,..., m} x {1,..., B} Z+ 

> Single-photon detection counts. 

• Be Z + > Number of times single-photon counting is repeated . 

Output: 

• {&ij} G (R + U 0) m > Estimate of {a t j} (8). 

1: procedure AMPLlTUDEESTlMATlON(m, N ijbj , B) 

2: for i,j G {1, • • •, m} x {1,..., m} do 

3: otij <— Mean( y/N nbl N ijbj JN ljbj N ilbl '■ b±,bj G {1,... ,-B}) 

4: end for 

5: return {ay,-} 

6: end procedure 


The probabilities Pij are estimated by counting single-photon detection events. 
Suppose Nij b . photons are transmitted from input port j to the detector at output 
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port i when N bj photons are incident and bj E {1,... ,B}. For large enough B, the 
transmission probability converges according to 

Pij mean j : bj E {1,..., £} j . (17) 

Likewise, the amplitudes {cep} are estimated by averaging the single-photon detection 
counts according to 


Oiij 


' Pnpjj 

PljPil 


<— mean 


Nijbj Nbj N bl 

Nb! Nbj Nljbj Nii bl 


: bi, bj E {1, • • •, B} 


= mean 


' NnfoNijbj 
NijbjNnb! 


: bi, bj E {1, • • •, B} 


(18) 


The estimate of cep relies on single-photon counts measured by impinging photons at 
the first input port repeatedly (repetition index b\ G {1,..., B}) and independently at 
the j’-th input port (with repetitions labelled by a different index bj G {1,..., B}). 

Henceforth, we represent our estimate of any parameter • by •. The estimate 
otij calculated using (18) is independent of N bj and thus resistant to variations in the 
incident-photon number N bj over different input modes j and different repetitions bj. 
Thus, our estimates {dp } are accurate in the realistic case of fluctuating light-source 
strength and coupling efficiencies. 

Finally, the standard deviations a (dp) of our estimates are calculated according to 


u(dp) <— std. dev. : 6 i» b j E {1,..., , (19) 

which converges for a large enough B. In line with standard nomenclature, we refer to 
these standard deviations as error bars. Algorithm 2 details the estimation of {dp} and 
error bars on the obtained estimates. 


3.3. Calibration to estimate mode-matching parameter 7 (Algorithm 3) 

In this subsection, we describe the procedure to calibrate our light sources for imperfect 
mode matching. The mode-matching parameter 7 is estimated using one- and two-photon 
interference on an arbitrary beam splitter. First, the reflectivity of the beam splitter is 
determined using single-photon counting [33]. Next, controllably delayed photon pairs 
are incident at the beam splitter inputs and coincidence counting is performed on the 
beam splitter output . We introduce a curve-fitting procedure to estimate the value of 7 
such that ( 12 ) best fits the measured coincidence counts. 

The beam-splitter reflectivity, which is denoted by cost?, is estimated as follows. A 
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which is in the form of (8) with «22 = cot 2 t). The value of «22 is estimated using single¬ 
photon counting as described in Algorithm 2. The estimated beam-splitter reflectivity 



The error bar on cos d is estimated by repeating the photon counting along the lines of 
Algorithm 2. 


Algorithm 3 Calibration Calculates the mode-matching parameter 7 of source- 
held using a beam splitter of known reflectivity. 

Input: 

• k,Q — {oj\, o> 2 ,... Uk-i, 0Jk} £ (R + ) fc > Frequencies at which / 1 , / 2 are given. 

• / 1 , fi '■ 12 —i > R + > Given spectral functions. 

• £,T — {ti, T 2 , ■ • •, 77 } G (R U 0)V Time delay values coincidence is measured at. 

• C cal : T —» R + > Measured coincidence curve. 

• i) G (—7T, 7r] > cos d is rehectivity of calibrating beam splitter. 

Output: 

• 7 G [0,1] > Estimate of mode-matching parameter of photon source. 

1 : procedure Calibration^, hi, / 1 , / 2 , £, T, C cal , i)) 

2: A G- {cost?, sin d, sin i ) 1 cost?} > Beamsplitter of rehectivity R (20) 

3: G- {0, 7t/2, 7t/2, 0} > Beamsplitter of rehectivity R (20) 

4: C(t, 7) = Coincidence ( 12 , / 1? / 2 , T, A, <f>, 7) > The quantities 12 , fi, / 2 , R cal are 

given. 7 is unknown. Coincidence(12, /1, /2, T, i? cal , 7) depends on 7 and r 
5 : return 7 Fit(C(t, 7), C cal (r), 1 /C cal (r), InitGuesses) > Least-squares curve 

IC ca i( t)_ C(t 'y)P 

htting to obtain the value of 7 that minimizes ——— gcal ^ - 1 —. The argument 

1/C cal (r) is the weight function [44] that accounts for experimental noise, which 
is assumed to be proportional to \JC(r). Ignore values of r at which C(t) = 1. 
Appendix B details the choice of initial guesses to the algorithm. 

6 : end procedure 

Next we estimate 7 using two-photon coincidence counting. Controllably delayed 
pairs of photons are incident at the two input ports of the beam splitter. Coincidence 
measurement is performed at the output ports for different values of time delay between 
the two photons. A curve-fitting algorithm is employed to find the best-ht value of 7, i.e., 
the value 7 that minimizes the squared sum of residues between the measured counts and 
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the coincidence counts expected from (12) for the beam splitter matrix (20). Algorithm 3 
details the calculations of 7, which is used to estimate {%} values accurately. 


S.f. Two-photon interference to estimate {dij} (Algorithms 4-6) 


In this subsection, we describe our procedure to estimate the arguments {dij} of the 
representative matrix U ( 8 ). Our procedure requires the measurement of coincidence 
counts for 2 (m — l ) 2 different choices of input and output ports. Of these measurements, 
(m — l ) 2 are used to estimate the absolute values { 1 d l3 |} of the arguments and the 
remaining (m — l ) 2 are used to estimate the signs {sgn#p}. 

The absolute values {1 d l3 |} are estimated as follows. Single-photon pairs are incident 
at input ports 1 and j and coincidence measurements are performed at output ports 
1 and i for i , j G {2, ... , m}. The state (3) of a photon pair is transformed under the 
action of the 2x2 submatrix 


Uiljl 


y/Kl 0 \ 

0 Hr J 


A o)('i 1 \ A/mT o \ o \ 

0 ^ 0 ^TjJ ( 0 


( 22 ) 


of U labelled by the rows 1 and i and columns 1 and j. The probability of detecting a 
coincidence at the output ports 1 ,i is 


Cnji^r) =K j K 1 \ j \ 1 V i V 1 HiHi (a 2 - + 1} J d(Jid(J2\fj(uii)fi(uJ2)\ 2 

+ 27 Qij [ d^d^ fj(ui) ffjM fiM cos (La 2 T-LaxT + Oij) , (23) 


which is obtained by setting i! = j' = 1 in ( 12 ). 

The measured coincidence counts are used to estimate the value of d , 3 ] as follows. 
The shape of the coincidence-versus -t curve (23) depends on the values of a l3 and %. 
The shape does not depend on the parameters K t , Ai, A*, pi, pj, ui, v v which lead to 
a constant multiplicative factor to the coincidence expression. Furthermore, the shape 
is unchanged under the transformation —» — % for d %3 G (— 7 r, 7 r] if the spcct-ral 

functions are identical. Hence, 1 6 i3 \ can be estimated using the shape of the coincidence 
function (23) and the values {dp} estimated using Algorithm 2. A curve-fitting algorithm 
estimates the value 0 t3 \ G [0, 7 r] that best fits the measured coincidence counts. The 
calculation of {l^pl} is detailed in Algorithm 4. 

Our procedure computes the signs by using an additional (m — l ) 2 coincidence 
measurements. First we arbitrarily set 622 as positive 


sgn 622 = 1 (24) 

because of the invariance|| of one- and two-photon statistics under complex conjugation 
U —y U* [33]. The signs of the remaining arguments {dij} are set using the coincidence 

|| Expectation values of Fock-state projection measurement with Fock-state inputs are unchanged 
under U —> U* if the spectral functions are equal /i(w) = Afw). Otherwise, the sign of —022 can be 
ascertained using the difference in the r > 0 and r < 0 coincidence counts in 02 , 2 , 1,1 (t). 
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Algorithm 4 Argument2Port: Calculates the unknown complex argument in 
the entries of a 2 x 2 transformation using a two-photon coincidence curve. 

Input: 

• k, ft — {o;i, o> 2 , • • • Wfc_i, ujk} G (R + ) fc > /i, f 2 are measured at frequencies Q. 

• fi, f 2 : —>■ R + > measured spectra. 

• £, T — {ri, t- 2 , • • •, Tf\ G (R U 0)V Time delay values coincidence is measured at. 

• C exp : T —> R + > Measured coincidence curve. 

• A 4— {oiij, otiji,ai>j, oii'j'} > Complex amplitudes of 2 x 2 submatrix of A (8). 

• 0f- G (—7r, 7r]} > Three complex arguments of submatrix. 

• 7 > Mode-matching parameter of photon source. 


Output: 

. \e. 


yi > Estimated magnitude of the unknown complex argument, 

procedure Argument2Port(/c, fl, fi, / 2 , £, T, C exp , A, 0, 7 ) 

<f> = {9ij, 9ijt, 9j/j, 9i'j'} > Set of three known phases and one unknown phase. 

C(t, 9 if) = Coincidence (hi, / x , f 2 , T, A, $, 7 ) 


return 9 i3 t— |LM(C(r, 0,j), C exp (r), 1/C exp (r))| 


> Use curve fitting to compute the value that minimizes 

5: end procedure 


E1t6T |Cexp(T)— C(r,7)| 2 

C*exp (t) 


counts between output ports when photon pairs are incident at input ports {j, j' j 

for a suitable choice of as we describe below. The coincidence probability at the 

output ports i, if is 


Ca'jjfr) =K i Ki>\i\i>n j n j 'V j v j > [ + a^alf) j dcoi dco 2 1 /y (coi) (co 2 ) | : 

+ 2'ya ij a ij 'a i > j a i ' j > [ du 1 du 2 f j (uj 1 )f f (uj2)f j (u 2 )f j '(uj 1 ) 


X COS (uj 2 T — U>iT + Pii'jj 1 ) 


(25) 


where 

flu'jj' — 1 9i'ji — 9ijf — 9i'j + 9ij\ e [ 0 , 7 r]. (26) 

Curve fitting is employed to estimate the value of fin> 33 ’ that best fits the measured 
coincidence counts. 

The estimated value of fiu'j-j’ is employed by Algorithm 5 to ascertain the sign of 
9ij. Algorithm 5 relies on the identity 


sgn 9ij = sgn (| p Wjj , - P iVjj ,\ - \p iVjj , - p? Vjj ,\) , 
and on known values of 


Aw = 1 9tf-9 if - 0,,j ± Iflyll,Atw e [0,7T 


(27) 


(28) 


to ascertain the sign of 9 VJ . If the sign of % is positive, then f u ' 3J ’ = and (27) 

returns a positive sgn 9 i3 . Otherwise, f u ' 33 ' = •■,, in which case (27) gives a negative 

sign. 
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Algorithm 5 SignCalc: Calculates the complex-phase sign of an element of the 
2x2 submatrix of an interferometer transformation matrix. 

Input: 



+ 

oT 

1 

|^4 

1 

^~4 

hi 

• 

> As defined in (27). 


• 0#^, 9ijt, 6i'j, \@ij\ 

> Equations (27-28). 

Output: 



• sgn % 

> Sign of 9 G (—7r, 7r] is defined in (14) 

1 

procedure SignCalc(/3, 9pj>, 6^, 9pj , % ) 


2 

(5 + \9i'ji — 9{ji — 9i>j + \9if\ 

> If 9 g h < 0, then /3 — /3~. 

3 

(j t 9}/ j ! f'Cj' ^i'j ffp’ 

> If 9 gh > 0, then fd = j3 + . 

4 

sgn 9 i:j G- sgn \/3 - /3~ \ - \/3 - fd + \ 


5 

return sgn^j 


6 

end procedure 



Algorithm 6 iteratively chooses indices i,i',j,f such that the signs of 9 t] p Opj, 9,^' 
have already been ascertained before ascertaining the sign of 9 l3 . In each iteration, the val¬ 
ues of are calculated by substituting 9 l3 |, |%'|, \0i r j\, dry |, sgn0 l3 >, sgn 9,p , sgn 9,pi. 

The algorithm estimates fja'jf by curve fitting measured coincidence counts to (25). 
Algorithm 5 is ascertains the sign of % using the estimates of ff n 'jp and One 

suitable ordering of indices ii'jj', which we depict in Figure 5, is 

• set if = 2, j' = 1 to determine sgn#j 2 for i G {3,... pm} (Figure 5b), 

• set if = 1 ,j' = 2 to determine sgn0 2 j for j G {3,..., m} (Figure 5c), 

• set if = 2 ,j' = 2 to determine sgn 9^ for (i,j) G {3,..., m} x {3,..., m} (Figure 5d). 

In summary, sgn6b,- is determined using the values of Pnpp , which are estimated by curve 

fitting, and of , which are computed using the signs and amplitudes of 0 l3 >, 9 t > 3 , 60/ • 
Algorithms 4-6 detail the step-by-step procedure to determine the absolute values and 
the signs of {%}. 

For certain interferometers U, the ordering of indices ii'jj' depicted in Figure 5 
can lead to instability in the characterization procedure. Appendix A elucidates on 
this instability and presents strategies to counter the instability. This completes our 
procedure to characterize the matrix A for representative matrix U = LAM. In the next 
subsection, we present a procedure to estimate the matrix that is most likely for the 
characterized matrix A. 

3.5. Maximum-likelihood estimation for finding unitary matrix 

At this stage, we have estimated the matrix A (8). The diagonal matrices L and M 
can be uniquely determined from A as follows. The representative matrix U = LAM is 
unitary so we have 


f/C/t = 1, 


( 29 ) 
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Figure 5: A depiction of the sign estimation procedure in Lines 9-22 of Algorithm 6. 
(a) The first row and first column arguments {@a}, {@ij} are zero, so their signs are 
arbitrarily set as positive. 9 22 is set as positive according to (24). (b) The sign of 
each second row argument 9 i2 is set using the known values |^ 22 1, \Qa\ and coincidence 
measurement for input ports 1,2 and output ports 2 ,i as in Line 15. (c) The sign of each 
second column argument 0 2 j is set using the known values \9 22 \, d‘ 2 j | and coincidence 
measurement for input ports 2, j and output ports 1,2 as in Line 16 of Algorithm 6. 
(d) The signs of each remaining argument 9 l3 is set using the known values \0 22 \, 0 l2 |, \9 2 j\ 
and coincidence measurement for input ports 2 ,j and output ports 2 ,i in Line 19 of 
Algorithm 6. 


which, upon substitution U = LAM , implies that 

LAMM* A 1 L* = 1 
=► AMM*A ] = L~ l L*~ l . 


Considering the first columns of the matrices (30) gives 




n\ 

0 


or 


A 


Similarly, using U'U = 1 we obtain 


( Li\ 

\Lm / 


\0 / 


A\ 

0 

Vo ) 



f 1\ 


n\ 

A f 

^2 

1 

0 


\\n/ 

hi 

w 


(30) 


(31) 


(32) 


( 33 ) 
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Algorithm 6 ArgumentCalc: Calculate {0 tJ } using two-photon coincidences 

Input: 

• k, Cl = {u)i,u> 2 , ■ ■ ■ u>k-i,uJk} G (R + ) fc > fi, f 2 are measured at frequencies Cl. 

• fi, f 2 '■ Cl R + > measured spectra. 

• £,T — {ri, T 2 ,..., 17 } G (R U 0)u> Time delay values coincidence is measured at. 

• C WjA T ) for (hi'Jif) e I 1 ’ 2 ) x { 1 ,... ,m} x { 1 , 2 } x { 1 ,.. .pm}, i ± iff ^ f 

> Measured coincidence at output ports if f when photons that have mutual 
delay r are incident at input ports i,j. 

• a : {2,, m} x {2,..., m} —» R + > Complex amplitudes ( 8 ). 

• 7 G [0,1] > Mode-matching parameter estimated using Algorithm 3. 

Output: 

• 9ij : { 1 ,..., m} x { 1 ,..., m} —> (— tt, i r] > Complex Arguments ( 8 ). 

1 : procedure ArgumentCalc(/c, Cl, fi, f 2 , £, T, CffCfr), a tJ , 7 ) 

2 : for i in { 1 ,..., m} do 

3: On, On, sgnflji, sgndij ^—0 > The first row, column are real valued. 

4: end for 

5: for (i, j ) in {2,..., m} x {2,..., m } do 

6 : A G- {1,1,1 ,oi g h}, {0,0,0} >2x2 matrix: rows 1 ,i, columns 1 ,j. 

7: \6ij | t— Argument2Port (C^jT, Cl, f u f 2 , T, A, d>, 7 ) 

8 : end for 

9: sgn 0 2 2 1 > The sign of 0 22 is positive by definition. 

10: for (ijifjjf) G {2}x{3,... ,m}x{2}x{3,... ,m}U{l}x{2}x{2}x{3,... ,m}U 

{2} x {3,..., m} x {1} x {2} do 
11: A G- {0,0,0}, G- {0,0,0} 

12 : /Vjj' <- Argument2Port (C-*C,(t),CI, / 1 , f 2 , T , A, $, 7 ) 

13: end for 

14: for i in {3,..., m } do 

15: 0i 2 g- 1 0i 2 1 Sign Calc (/?i 22 * ? 0 ? ^ 22 ,0, | ^ 2 1,)); 

16: 0 2 i |02i|SlGNCALC(/32il2, 0, 0 2 2, 0, \0 2 i \,)); 

17: end for 

18: for {i, j ) in {3,..., m} x {3,..., m} do 

19: dij g- 1 0ij | Sign Calc (p Wjjl , 0 22 , 0 i2 , 0 2j , | % |) 

20 : end for 

21 : return {%} 

22 : end procedure 


Equations (32) and (33) are systems of linear equations that can be solved for L and M 
respectively using standard methods [45]. The solutions L and M of the linear systems 
and the characterized matrix A give us the representative matrix U = LAM. 

The experimentally determined A is different from the actual A because of random 
and systematic error in the experiment, where we denote the experimentally determined 
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Algorithm 7 MaxLikelyUnitary: Calculates unitary matrix that has maximum 
likelihood of generating estimated {ay,}, {%} 

Input: 

• a : {1,..., m} x {1,..., m } —> R + U 0 > Estimated amplitudes of A (8). 

• 9 : {1,..., m} x {1,, m} —> (— tt, 7r] > Estimated arguments of A (8). 

Output: 

• W G SU(n) > Unitary matrix with maximum likelihood of generating A. 

1: procedure MAXLIKELYUNITARY^-, %) 

2: Aj i — 1 

3: {/ij : i e {1,..., m}} solution of system (32) of linear equations. 

4: {Aj : i e {2,..., m}} solution of system (33) of linear equations. 

• r >. Up t \jCXjjC ' 1 fij 

6: W <— ( 'uU ^ 2 U > Assumption: U t j — Ujj is an iid Gaussian random variable 

with zero mean for all i, j G {1,..., m}. 

7: end procedure 


values of interferometer parameter • by •. Similarly, the L and M matrices obtained 
by solving Eqs. (32) and (33) for A (rather than A) differ from the actual L and M 
respectively. The estimated U = LAM is thus a non-unitary matrix and is not equal to 
U in general. Furthermore, U is a random matrix, which depends on the random errors 
in the one- and two-photon experimental data. 

We employ maximum-likelihood estimation to calculate the unitary matrix W that 
best fits the collected data. First, bootstrapping techniques are used to estimate the 
probability-density function (pdf) of the entries of the random matrix U [46,47]. Next, 
standard methods in maximum-likelihood estimation [48] are employed to find the unitary 
matrix W. Maximum-likelihood estimation simplifies under the assumption that the 
error on U is a Gaussian random matrix ensemble, i.e, that the matrix entries jCp j 
are complex independent and identically distributed (iid) Gaussian random variables 
centred at the correct matrix entries. In this case, the most likely unitary matrix W is 
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the one that minimizes the Frobenius distance^ from U [49]. The unitary matrix 



minimizes the Frobenius-norm distance from U [50]. Thus, if the random errors {U^ — U,^} 
in the matrix elements are iid Gaussian random variables with mean zero, then W is the 
best-fit unitary matrix. Figure 6 is a depiction of the actual, the estimated and the most 
likely transformation matrices. Algorithm 7 computes W. 

This completes our procedure to estimate the most-likely unitary matrix W that 
represents the linear optical interferometer. In the next section, we present a procedure 
to estimate the error bars on the entries of the estimated representative matrix W 
accurately. 



Figure 6: A depiction of the error in reconstruction of the interferometer matrix U. The 
matrix U represents the unitary transformation effected by the interferometer. U is 
the complex-valued transformation matrix returned by the reconstruction procedure. 
Algorithm 7 returns W, which represents the unitary matrix that is most likely to have 
generated the data collected in the characterization experiment. 


The Frobenius norm of a matrix A mxm is defined as 




\ 


Ei 


The Frobenius-norm distance between matrices U and V is defined as 


(34) 


dist(U,V) = \\U-V\\ F = 


i E i^- - A, 

\ i,j=l 


2 

I 5 


( 35 ) 


and is a symmetric, positive-definite and subadditive distance function on the set of matrices. 
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In this section, we present a procedure to estimate the error bars on the matrix entries 
{Wij} of the characterized representative matrix W. The entries { W t] } computed by 
Algorithms 1-7 are random variables because of random error in experiments. Obtaining 
accurate error bars on these random variables is important for using characterized linear 
optical interferometers in quantum computation and communication. Current procedures 
compute error bars under the assumption that Poissonian shot noise is the only source 
of error in experiment [21,23]. 

We choose to employ bootstrapping on the data determine error bars [46, 47,51-53]. 
Monte-Carlo simulation is widely used but this technique is not applicable here because 
the Poissonian shot noise assumption is not reliable given the presence of other sources 
of error some of which are not understood. Bootstrapping is preferred because the 
nature of the error need not be characterized and instead relies on random sampling with 
replacement from the measured data. Bootstrapping can be employed toyield estimators 
such as bias, variance and error bars. 

Algorithm 8 calculates the error bars a(W t j) using estimates of the {Wy} pdf’s, which 
are obtained using bootstrapping as follows. The algorithm simulates N characterization 
experiments using the one- and two-photon data, i.e., the inputs to Algorithms 1-7. 
In each of the N rounds, the one- and two-photon data are randomly sampled with 
replacement (resampled) to generate simulated data. The data thus simulated are given 
as inputs to Algorithms 1-7, which return the simulated representative matrices 

{W ,b : b e {1,...,A}, N e Z + } . (37) 

The pdf’s of the simulated-matrix entries {W' b : b G {1,... ,1V}} converge to the pdf’s 
of the respective elements {} for large enough N [54,55]. 

The simulated data are obtained in each round by resampling from the one- and 
two-photon experimental data as follows. Single-photon detection counts are simulated 
by resampling from the set {N^ : bj e {1,..., B}} of experimental detection counts 
(Line 17 of Algorithm 8). Two-photon coincidence counts are simulated by shuffling 
residuals obtained on curve-fitting experimental data. Specifically, the algorithm (Line 12) 
resamples from the set 


Mr) = - Cn'jji (t) : t G T} (38) 

of residuals obtained by fitting experimentally measured coincidence counts to function 
Cu'jj>( T ) (12). The resampled residuals are added to the fitted curve to generate the 
simulated data (Line 14) + . Algorithms 1-7 are used to obtain the simulated elements 

+ The pdf of the residuals is different for different values of r. We assume that the pdf’s for different 
r are of the same functional form, albeit with different widths. The distribution of the residuals 
for different values of r are determined using standard methods for non-parametric estimation of 
residual distribution [56,57]. Algorithm 8 normalizes the residuals before resampling from the residual 
distribution. 



Accurate and precise characterization of linear optical interferometers 


22 


Wij of the representative matrix. Finally, the error bars on the { W t] } are estimated by 
the standard deviation of the pdf of the elements. 

This completes the characterization of representative matrix W and the error bars 
on its elements. The next section details a procedure for the scattershot characterization 
of the interferometer to reduce the experimental time required for characterizing a given 
interferometer. 

5. Scattershot characterization for reduction in experimental time 

In this section, we present a scattershot-based characterization approach to effect a 
reduction in the characterization time [58,59]. Our scattershot approach reduces the 
time required to characterize an m-niode interferometer from O (m 4 ) to O (m 2 ) with 
constant error in the interferometer-matrix entries. 

The straightforward approach of characterization involves coupling and decoupling 
light sources successively for each one- and two-photon measurement. In contrast, the 
scattershot characterization relics on coupling heralded nondetcrministic single-photon 
sources to each of the input ports of the interferometer and detectors to each of the 
output ports. Controllable time delays are introduced at two input ports, which are 
labelled as the first and second ports. All sources and detectors are switched on and the 
controllable time-delay values are changed first for the first port and then for the second 
port. 

Single-photon data are collected by selecting the events in which exactly one of 
the heralding detectors and exactly one of the output detectors register a photon 
simultaneously. Two-phot.on coincidence events at the outputs are counted when 
two heralding detectors register photons. The controllable time delays introduced 
at the first and second input ports ensure that each of the 2(m — l) 2 coincidence 
measurements is performed. Note that our characterization procedure (Algorithms 1-8) 
yields accurate estimates of interferometer parameters even when photon sources with 
different spectral functions are used. In summary, the required characterization data 
are collected by selectively recording one- and two-phot.on events. The setup for the 
scattershot characterization of an interferometer is depicted in Figure 7. 

Now we quantify the experimental time required in the characterization of a linear 
optical interferometer. Our characterization procedure requires Bm 2 single-photon 
counting measurements and 2(m — l) 2 coincidence-counting measurements to characterize 
an m-mode interferometer. We estimate the time required for each of these measurements 
such that random errors in the {0,-j} estimates remain unchanged with increasing 

m. To ensure constant error in the {ay,}, {@ij} estimates, we require that the number 
of one- and two-phot.on detection counts remain unchanged with increasing m. The 
probability of photon detection at. the output decreases with increasing m because of 
the concomitant decrease in the transmission amplitudes {«p}- 

The amplitudes {a^} drop as O (1 /\Jm) because of the unitarity of U [60]. Hence, 
one- and t.wo-phot.on transmission probabilities (10,12) decrease as 1/m and 1 /m 2 , 
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Algorithm 8 BOOTSTRAP: Estimate error bars on {Wij}. 

Input: 

• k, fi, /i, f 2 : fi —>• R + > Spectral functions: same as Algorithm 3. 

• £, T — {ti, t 2 , ..., Te} G (R U 0) £ > Time delay values. 

• m, C cal (r), C t yh,(r) for r G T and (■) > Same as Algorithms 3 and 6 . 

• B, Nijbj : {1,..., m } x {1,..., m} x {1,..., B} —» Z + as in Algorithm 2. 

• A > Number of bootstrapping samples. 

Output: 

• er (re Wy), cr(im W l3 ) : {1,..., m } x {1,..., m} —» R + > Error in W elements. 
1 : procedure Bootstrap^, fi, fi, f 2 ,£, T, 7 , C'^h,(r), %, B, N ijb , A) 

2: A G- {cos'd, sin d, sind, cos'd}, G- (0, 7 t/ 2 , 7 t/ 2 , 0} 

3: Residuals cal (r) G- C cal (r) - Coincidence^, Q, f u f 2 , £, T, A, 0, 7 ) 

4: NormalResiduals cal (r) G- Res1 ^^ (A > Assumption: Resicluals cal (r) pdf width 

oc C m (r). 

5: for (i,i',j,f) G {1,2} x {1,... ,m} x {1,2} x {1,.. .,m}, i ± i',j ± f do 

6 : A 4 {QTy, Oii'j) Olij ', O 7 }, 4 {d^y, Oi'j, Qij /, d^'} 

7: C^ y g- Coincidence^, fi, / 1 , / 2 , £, T, A, 0, 7 ) 

8 : Residuals*/^ (r) G- C£? jf (r) - Cf^y (r) 

9: NormalResiduals n'jj'(r) G- 

jj f ' ' 

10 : end for 

11 : for n = 1 to A do 

12 : ShuffledNormalResiduals cal (r) G- Resample |T| residuals from 

NormalResiduals cal (r) 

13: ShuffledResidual cal (r) G- C fit (r) x ShuffledNormalResiduals cal (r) 

14: C n {r) = C fit (r) + ShuffledResidual(r) 

15: 7 n G- Calibration^, fi, / 1 , / 2 , f, T, C b , 'd) 

16: for G {1,2} x {l,...,m} x {1,2} x {l,...,m}, i ^ / do 

17: -(—Mean y/An;,, A,y, /.Yy, Aab} from B values each of bi,bj drawn 

with replacement from {1,..., B} 

18: ShuffledNormalResidualsyy/ (r) G- \T\ entries in Nor malResiduals j jy/ (r) 

19: ShuffledResidualjj/jj/ (r) G- Cff, r] , (r) x ShuffledNormalResidualSjjy/(r) 

20 : C^ijji (t ) = (^ 77 /(t) + ShuffledResidualjj/jj^r) 

21 : end for 

22 : {d£} = ArgumentCalc(A;, fi, / 1 , / 2 , £, T, Qy, (t), cy, 7 ) 

23: {W^} = MAXLIKELYUNITARY({«”■}, {d-}}) 

24: end for 

25: for (i,j) in {1,..., m} x {1,..., m} do 

26: a(reWij) = std. dev. ({r eWij }); cr(inilly) = std. dev. ({imWjj}) 

27: end for 

28: end procedure 
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Figure 7: Schematic diagram of the procedure for scattershot characterization of a 
five-mode interferometer U. Heralded single-photon sources are coupled to the inputs of 
the interferometer and controllable time delays are introduced at the first two ports. All 
sources and detectors are switched on and the controllable time delay values are changed 
for the first port and then for the second port. The required characterization data are 
collected by selectively recording one- and two-photon events. 


respectively. More photons need to be incident at the interferometer input ports to 
offset this decrease in transmission probabilities. Therefore, maintaining a constant 
standard deviation in the {cpj} and {%} measurements requires O (m) and O (m 2 ) 
scaling respectively in the number of incident photons, which amounts to an overall 
O (m 4 ) scaling in the experimental time requirement. Scattershot characterization allows 
(m — l) 2 different sets of the one- and two-photon data to be collected in parallel thereby 
reducing the time required to characterize the interferometer by a factor of (m — l) 2 . 
The overall time required for the characterization decreases from O (m 4 ) to O (m 2 ) if 
the scattershot approach is employed. 

Our analysis of scattershot characterization assumes that the coupling losses are 
small and that weak single-photon sources are used, i.e., that the probability of multi¬ 
photon emissions from the heralded sources is small as compared to single-photon emission 
probabilities. These assumptions are expected to hold for on-chip implementations of 
linear optics that have integrated single-photon sources and detectors. 

Light sources used at each input port in our scattershot-based characterization 
procedure differ spectrally in generally. Our characterization procedure is accurate 
despite this difference because we measure source-field spectra and using these data in 
the curve-fitting procedure. 

We have developed the scattershot approach which has advantages and disadvantages 
but on balance is a superior experimental approach to consecutive measurement. The 
advantage is that the time requirement for characterization if reduced by a factor that 
scales as O (m 2 ). The disadvantage is the overhead of requiring one source at each input 
port and one detector at each output port. The disadvantage is not daunting because 
these requirements are commensurate with other active investigations of QIP such as 




























Accurate and precise characterization of linear optical interferometers 


25 


LOQC and scattershot BosonSampling. In fact, state of the art implementations [59] 
meet our increased requirements for scattershot characterization. 

6. Summary of procedure and discussions 

In this section, we summarize our characterization procedure for a less formally-oriented 
audience. We describe the processing of the collected experimental data by the various 
algorithms presented in Section 3. We compare our procedure with the existing procedure 
for the characterization of linear optics using one- and two-phot.ons [33]. We provide 
numerical evidence that our characterization procedure promises enhanced accuracy and 
precision even in the presence of shot noise and mode mismatch. 

The experimental data required by our procedure to characterize an m-mode 
interferometer includes the following one- and two-photon measurements. The number 
Nij b - (13) of single-photon detection events is counted at the j'-tli output port when single 
photons are incident at the i-th input port. This single-photon counting is repeated 
B times for each of the input ports and output ports, where B is chosen such that 
the cumulants of the set {N^. : bj G {1, ..., B}} converge. The single-photon counts 
{TVjjfe,} are received by Algorithm 2, which returns the {ay,} (8) estimates using Eq. (18). 

The spectral function fj(co) (2) of the light incident at each input port j is measured. 
This function is used by Algorithm 1 to calculate the expected two-photon coincidence 
curves using Eq. 12. Fitting experimental data to these coincidence curves yields an 
accurate estimate of the mode-matching parameter during calibration and the arguments 
{Oij} in the argument-estimation procedure. Thus, the spectral function fjipj) serves as 
an input to the algorithms for the estimation of the mode-matching parameter and of 
the arguments {%} (Algorithms 3-6). 

The mode-matching parameter 7 is estimated by performing coincidence 
measurement on a beam splitter that is separate from the interferometer but is constructed 
using the same material. First, we use single-photon data to estimate the reflectivity 
cost? of the beam splitter according to Eq. (21). Imperfect mode-matching changes the 
shape of the coincidence curve, and we find 7 by comparing the shapes of (i) the curve 
expected for reflectivity cost? and (ii) the curve obtained experimentally. The estimated 
beam splitter reflectivity, the measured spectra and the coincidence counts are received 
as inputs by Algorithm 3, which returns an estimate of 7. 

Algorithm 6 uses two-photon coincidence counts to estimate the arguments {%}. 
Coincidence counts are measured for the input ports j,j' and output ports i,i ' for the 
2 (m — l) 2 sets 

e I 1 , 2 } x {1,... ,m} x {1,2} x {1,... ,m}, i ± 1 , j ± j’ (39) 

of input and output ports. In other words, coincidence counts are measured for different 
choices of two input ports and two output ports, such that each of the choices includes 
(i) either the first or the second input ports and (ii) either the first or the second output 
port. Algorithm 6 receives as input the measured spectra, the {ay,} values estimated 
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by Algorithm 2, the 7 value estimated by Algorithm 3 and the two-photon coincidence 
data for the choice (39) of input ports. The algorithm returns the {%} estimates. 
The computed estimates of {a t j} and of {0,j} yield the representative unitary matrix 
W (36) that has maximum likelihood of describing the characterized interferometer 
(Algorithm 7). This completes a summary of our procedure for characterization of the 
interferometer. 

Algorithm 8 employs bootstrapping to find the error bars on the elements { W tJ } of the 
characterized unitary matrix. The bootstrapping procedure uses the experimental data 
that is received by Algorithms 1-7 and repeatedly simulates experiments by resampling 
from the experimental data. The number N of repetitions is chosen such that the pdf’s 
of the {Wjj} elements over many rounds of simulation converge. The error bars on 
the {Wij} elements are computed based on the estimated pdf’s of the elements. Our 
procedure thus enables the estimation of meaningful error bars on the characterized 
unitary matrix. 

Bootstrapping is employed to test the goodness of fit between the experimental 
curve and expected curves [61]. Experiments [11,36] can employ bootstrapping instead 
of the incorrect x 2 -conhdence measure to test if the data are consistent with quantum 
predictions or with the classical theory. 

Finally, we recommend a scattershot approach for reducing the experimental time 
required to characterize interferometers. The approach involves coupling heralded 
nondeterministic single-photon sources at each of the input ports and single-photon 
detectors at each of the output ports. All the sources and the detectors are switched 
on in parallel. Single-photon counts are recorded selectively as two-photon coincidences 
between the heralding detectors and the output detectors, while two-photon events 
are recorded when two heralding detectors and two output detectors record photons. 
Controllable time delays are introduced at the first and second input ports so coincidences 
between each of the 2(m — l) 2 choices (39) of input and output ports are recorded. The 
scattershot approach reduces the experimental time required to characterize an m-mode 
interferometer from O (m 4 ) to O (m 2 ). 

Now we compare and contrast our procedure with the Laing-O’Brien procedure [33]. 
Our procedure is inspired by the Laing-O’Brien procedure in that it employs (i) a 
‘real-bordered’ parameterization (8) of the representative matrix and modelling of linear 
losses at the interferometer ports, (ii) a ratio of single-photon data to estimate the 
complex amplitudes of the matrix elements and (iii) an iterative procedure that uses 
two-photon data to estimate the amplitudes of the complex arguments and to estimate 
the signs of the complex arguments. 

Our procedure differs from the Laing-O’Brien [33] procedure in that we use averaged 
value (18) of the ratio of single-photon detections over many runs rather than the ratio of 
averaged values. This difference ensures accuracy of our procedure even under fluctuation 
in the number of incoming photons. Such fluctuations might arise from fluctuations in 
pump strength of the single photon source or in the strength of coupling between photon 
source and interferometer. 
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Figure 8 : The fitting of coincidence data to curves obtained from spectral functions 
using ( 12 ) and to Gaussian functions. Coincidence counts are simulated using 
experimentally measured spectra. 


Another advance in our method is the curve-fitting procedure for estimating complex 
arguments of interferometer matrices. The Laing-O’Brien procedure requires coincidence- 
curve visibilities to estimate complex arguments a^. Whereas the Laing-O’Brien 
procedure recommends coincidence probabilities be measured at zero time delay and 
also at time delays large as compared to the temporal spread of the wave-packet, in 
practice, current implementations determine the visibilities by fitting experimental 
data to Gaussian curves [35,63-67]. These implementations are flawed because source 
spectra differ from Gaussian in general. Our procedure is accurate because the data 
are fit to curves computed from spectral functions, rather than fitting to Gaussians. 
Figure 8 illustrates the distinction between fitting experimental coincidence counts to 
the coincidence function (12) simulated using spectra and fitting to Gaussian functions. 
Figure 9 demonstrates the increase in accuracy and precision of characterization by using 
the correct curve-fitting function. 

We introduce the calibration subroutine, which relies on the estimation of the mode 
mismatch in the source field. Spatial and polarization mode mismatch is not an issue 
of major concern in waveguide-based interferometers, which typically operate in the 
single-photon regime. In these interferometers, the calibration step of our procedure can 
be neglected without decreasing accuracy. The mode-mismatch parameter 7 , which is 
an input of the curve-fitting procedure, is set to unity. 

In the context of bulk-optics, our calibration step ensures accuracy and precision if 
(i) 7 is identified as the maximum-possible source overlap in the spatial and polarization 
degrees of freedom and (ii) the experimentalist adjusts the setup to maximize coincidence 
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Effect of fitting-curve choice on accuracy of characterization 



Figure 9: A plot showing the effect of fitting-curve choice on the accuracy and precision 
of the characterization procedure. The two curves depict the mean error for the two 
different choices of fitting curves, where the error is the trace distance between the 
expected and the actual unitary transformations and the mean is over 1000 simulated 
characterization experiments. One- and two-photon interference data was simulated 
for a five-channel interferometer using experimentally measured spectra and simulated 
Poissonian shot-noise. Characterization was performed by fitting coincidence curves to 
Gaussians (red curve) and to correct curves according to our procedure (blue curve). 
MATLAB code for the simulations depicted in this figure is available on GitHub [62] 


visibility for the calibrating beam splitter and for each choice of interferometer inputs 
ports. Such an adjustment will ensure that the source overlap acquires its maximum- 
possible value 7 in each of the coincidence-curve measurements. This maximum value is 
a property of the sources used and is independent of source alignment and focus so is 
expected to remain unchanged between different confidence measurements. Figure 10 
demonstrates the increase in accuracy and precision of characterization by using the 
calibration procedure 

Other advances made in our characterization procedure over existing procedures 
include (i) a maximum likelihood estimation approach to determine the unitary matrix 
that best fits the data (ii) a bootstrapping based procedure to obtain meaningful 
estimates of precision and (iii) a scattershot-based procedure to improve the experimental 
requirements of characterization. 
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Effect of calibration on accuracy of characterization 



Figure 10: A plot showing the effect of calibration on the accuracy and precision of the 
characterization procedure. The two curves depict the mean error for characterization 
with (blue curve) and without calibration (red curve), where the error is the trace distance 
between the expected and the actual unitary transformations and the mean is over 1000 
simulated characterization experiment. The simulations comprised generating one- and 
two-photon interference data based on experimentally measured spectral functions and 
performing characterization by our procedure, matlab code for the simulations depicted 
in this figure is available on GitHub [62] 


7. Conclusion 

In conclusion, we devise a one- and two-photon interference procedure to characterize 
any linear optical interferometer accurately and precisely. Our procedure provides an 
algorithmic method for recording experimental data and computing the representative 
transformation matrix with known error. 

The procedure accounts for systematic errors due to spatiotemporal mode mismatch 
in the source field by means of a calibration step and corrects these errors using an 
estimate of the mode-matching parameter. We measure the spectral function of the 
incoming light to achieve good fitting between the expected and measured coincidence 
counts, thereby achieving high precision in characterized matrix elements. We introduce 
a scattershot approach to effect a reduction in the experimental requirement for the 
characterization of interferometer. The error bars on the characterized parameters are 
estimated using bootstrapping statistics. 

Bootstrapping computes accurate error bars even when the form of experimental 
error is unknown and is, thus, advantageous over the Monte Carlo method. Hence, 
our bootstrapping-based procedure for estimating error bars can replace the Monte 
Carlo method used in existing linear-optics characterization procedures. We open the 
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possibility of applying bootstrapping statistics for the accurate estimation of error bars 
in photonic state and process tomography. 

Acknowledgments 

We thank Matthew A. Broome, Jens Eisert, Dirk R. Englund, Sandeep K. Goyal, Hubert 
de Guise and Michael J. Steel for useful comments. ID and BCS acknowledge AITF, 
China Thousand Talent Plan and NSERC and USARO for financial support. AK 
acknowledges CryptoWorks21 and NSERC for financial support. 



Accurate and precise characterization of linear optical interferometers 

Appendix A. Removal of instability in characterization procedure 


31 


In this section, we describe an instability in our characterization procedure, which can 
yield large error in the { W tJ } output for small error in the experimental data (r) in 
case of certain interferometers W. We present a strategy to circumvent this instability 
by means of collecting and processing additional experimental data. 

The instability in the characterization procedure arises because of an instability 
in estimation of {sgn dy\ (Algorithm 5). Small error in the measured coincidence 
counts can lead to the wrong inference of sgn 6^-, which can lead to a large error 
\\W — U\\ in the characterized matrix W. Recall that Algorithm 5 uses the identity 
sgn % = sgn - Pti'jyl - | fin'jj' ~ fih< 3 y\) (27) to determine the sign of the 

arguments, where fi^jy = |Oyy — diy — d^j ± |% 11 and the values of fiu'jy , Qyy , d^j , diy , |% \ 
are estimated by curve fitting. 

Random and systematic error in measured coincidence counts can lead to estimate 
of variables finyy , d,jy , diy , dy , 1 dy differing from their actual values. The estimation of 
sgn dij is unstable if the dyy — d t y — dy term (27) is close to 0 or tt because, in this case, 
a small error in the fiuyy estimate can lead to an incorrect sgn dy estimate. In other 
words, the sign estimates are unstable if the values of 


djj/jj* min dij' dj/j , r idi'j' dij> dj/jf , 


(A.l) 


are small compared to the error in our Aipy, d^y, d 7 >j, d t j/, \dij\ estimates. 
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Figure Al: An illustration of the instability in the d^ characterization procedure for 
an interferometer with m — 5 modes, (a) If the value of |do 2 — $52 — ^ 24 1 is close to 
0 or 7T, then small error in (r) can lead to an error in the estimation of sgn(0 54 ). 
(b) The instability in the #54 can be removed by collecting two-photon coincidence data 
for output ports 2, 5 and input ports 3,4 and using the values of 6 * 23 , d 53 , d 2 4 instead of 
& 22 , # 52 , d 2 4 values. 
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We mitigate the sign-inference instability by making two modifications to our 
characterization procedure; the first- modification removes instability from the sign- 
inference of the second row and second column elements whereas the second modification 
prevents incorrect inference of the remaining signs. The inference of {sgn 6 ** 2 }, {sgn 02 j} 
(Lines 14-17, Figures 5b, 5c) is unstable if 

0&j2 = min( 6 » 22 ,7r - 0 22 ) (A.2) 

is small as compared to the error in the /3 i2 j 2 , 0 22 , 0 2 j , 9 i2 , \6ij\ estimates. Hence, we relabel 
the interferometer ports such that 6 22 is as far away from 0 and n as possible. Specifically, 
after the amplitudes of the phases have been estimated (Line 8 of Algorithm 6 ), we 
choose i,j for which 1 9^ — 7t/ 2| is minimum, and we swap the labels of input ports 
2 , j and output ports 2 ,i. We measure two-photon coincidence counts based on this 
new labelling and process it using Algorithm 6 . The instability in the procedure for 
estimation of the {&&}, {0 2 j} signs is removed as a result of the relabelling. 

The second modification is aimed at removing the instability in the remaining 
signs. The procedure estimates the remaining signs by using {C(yL, (r)} values for 
i' = j' = 2. The estimation of 9ij is unstable if 9 x ff ]2 is small as compared to the error in 
the /3 i2 j 2 ,9 22 ,9 2 j,9 i2 , 9 l} estimates. We make a heuristic choice of a threshold angle 9 T 
t-hat accounts for the error in these variables, and we reject any sgn 0 tJ inferred using 
9j 2 j 2 < 9 T . Additional t-wo-photon coincidence counting is performed and employed to 
estimate these values of 0 tJ , as detailed in the following lines that can be added to the 
algorithm to remove the instability 
17 + 1 : for (i, j) in {3,..., m} x {3,..., m} do 
2: if 0 r l2j2 < 9 t then 

3: Choose i' ^ l,i and j' ^ 1 , j such that O'j' — %' — 9 y f is closest to ir/2. 

4: C^jfr) <7- Coincidence counts for input ports j,j' and output ports i,i'- 

5: A {0, 0, 0}, <f> <- {0, 0, 0} 

6 : Ai'ii'<- Argument2Port(C^v(t), ft, /i, / 2 , T, A, 4>, 7 ) 

7: 9ij <- % Sign Calc (0 ^,, 9 vy , 9 iy , 9 Vj , | % |) 

8: end if 

9: end for 

Figure A1 illustrates the rejection of those i,j choices for which 9f 2 y 2 < 9 T and the use 
(■ r),j' 7 ^ 2 counts to obtain a correct estimate of sgn 9 t j. We thus remove the 
instability in the estimation of {9^} and in the estimation of the representative matrix 
W. 

Appendix B. Curve-fitting subroutine 

Our characterization procedure employs curve fitting in Algorithm 3 to estimate the 
mode-matching parameter 7 and in Algorithms 4-6 to estimate {%} values. The curve¬ 
fitting procedure determines those values of unknown parameters that maximize the 
fitting between experimental and expected coincidence data. In this section, we describe 
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Figure Bl: Simulated coincidence counts for output ports i,i' and input ports j,f of 
interferometer with an = any = \/3/4 and a l% ' = cqy = 1/4 and for different values of 
Pn'jj’- The value of Pn'jj' in each respective figure is (a) 7 r, (b) 0, (c) tt/3 and (d) 2tt/3. 
The coincidence counts corresponding to r = 0 and r —)■ oo are marked on each plot by 
C' exp ( 0 ) and C exp (oo) respectively. 


the inputs and outputs of the curve-fitting subroutines. We present heuristics to compute 
good initial guesses of the fitted parameters. 

The curve-fitting subroutine receives as input (i) the choice of parameters to be fitted; 
(ii) the coincidence counts {C/ e /h,(r)}; (iii) an objective function, which characterizes 
the least-square error between expected and experimental counts; and (iv) the initial 
guesses for each of the fitted parameters. The output of the curve-fitting subroutine is 
the set of parameter values that optimize the objective function. 

The first input to the subroutine is the choice of the parameters to be fit. The 
curve-fitting subroutine fits three parameters. One of these three (namely the mode¬ 
matching parameter 7 in Algorithm 3 or the 1 9if or value in Algorithm 6 ) is related 
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to the shape of the curve, whereas the other two are related to the ordinate scaling and 
the abscissa shift of the curve respectively. The ordinate scaling factor comprises the 
unknown losses { K t , Uj }, transmission factors { X ,, pj } and the incident photon-pair count. 
The horizontal shift factor accounts for the unknown zero of the time delay between the 
incident photons. The algorithm returns the values of the shape parameter, the abscissa 
shift and the ordinate scaling that best fit the given coincidence curves. 

The objective function quantifies the goodness of fit between the experimental data 
and the parameterized curve. We use a weighted sum 

5>(t)|C“p(t)-C'(t)| 2 (B.l) 

rST 

of squares between the experimental data and the fitted curve as the objective function [44] 
for weighs w(r). We assume that the pdfs of the residues are proportional to ^/C exp (r) 
and we assign the weights 


w(t) 


1/C' exp (r) if C exp (Y) fz o 
1 if C ex p (r) = 0 


(B-2) 


to the squared sum of residues. In case the pdf’s of the residuals for different values of r 
is not known, standard methods for non-parametric estimation of residual distribution 
can be employed to estimate the pdf’s [56,57]. Thus, the curve fitting algorithm returns 
those values of the fitting parameters that that minimize weighted sum of squared 
residues between experimental and fitted data. 

The curve-fitting procedure optimizes the fitness function over the domain of the 
fitting parameter values. Like other optimization procedures, the convergence of curve 
fitting is sensitive to the initial guesses of the fitting parameters. The following heuristics 
give good guesses for the three fitting parameters. We guess the ordinate scaling as the 
ratio 

oo) 

Ca'jj' (oo) 

of the experimental coincidence counts 


(B.3) 



c «p ( oo) « gyM + 


(B.4) 


to the coincidence probability (oo) for large (compared to the temporal length 
of the photon) time-delay values. The 7 value is guessed for Algorithm 3 as the ratio 
of the visibility of the experimental curve to the expected visibility in the curve. The 
initial guesses for d = |dqj and d = f3u'jj> are based on the known estimate of 7 and the 

visibility 

T j. 2y cos 2 d sin 2 d . _ 

^ = 4 0 , • 4 Q • (B-5) 

cos 4 it + sm it 

of the curve. As there are four kinds of curves (see Figure Bl) possible for different 
values of the shape parameter (7, |%|, Ai'j'), another approach is to perform curve fitting 
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four times, each time with a value from the set 7 t/4, 37t/4, 5n/4, 7n/4 of initial guesses 
and choose the fitted parameters that optimize the objective function. Finally, the 
initial value of the abscissa shift parameter is guessed such that the global maxima or 
minima (whichever is further from the mean of the coincidence-count values over r) of 
the coincidence curve is at zero time delay. 

In summary, the curve fitting procedure uses the measured coincidence counts, the 
objective function and the initial guesses to compute the best fit parameters. This 
completes our description of the curve-fitting procedure and of heuristics that can be 
employed to computed the initial guesses for the fitted parameters. 
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